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In this paper we consider a microscopic model of a simple ecosystem. The basic in- 
gredients of this model are individuals, and both the phenotypic and genotypic levels 
are taken in account. The model is based on a long range cellular automaton (CA); 
introducing simple interactions between the individuals, we get some of the complex col- 
lective behaviors observed in a real ecosystem. Since our fitness function is smooth, the 
model does not exhibit the error threshold transition; on the other hand the size of total 
population is not kept constant, and the mutational meltdown transition is present. We 
study the effects of competition between genetically similar individuals and how it can 
lead to species formation. This speciation transition does not depend on the mutation 
rate. We present also an analytical approximation of the model. 

Keywords: Speciation models; Darwinian Theory; Population Dynamics; Eigen Model; 
Mutational Meltdown. 

1. Introduction 

Real ecosystems present a complex behavior. Many of their peculiar features 
are considered in classical population dynamics models, where the dynamical 
variables are the number of individuals of different populations £1 In this paper a 
different point of view is considered, we propose a microscopic model of an evolv- 
ing (in Darwinian sense) ecosystem, where the individuals are represented 
their genotypes. Our model is related to the Eigen model for quasispeciesp 
although we consider a different fitness landscape and the presence of interac- 
tions among individuals. 
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Considering simple ecological interactions (competition, predation, cooper- 
ation), we are able to obtain a complex collective behavior. The aim of this 
work is to have a simple predictive model, that can reproduce some of the basic 
features present in an real ecosystem, such as: 

• Evolution. The system has to be able to create diversity (mutations). 
Darwinian selection acts on this diversity. 

• Population dynamics. The model should reproduce the typical phenomenol- 
ogy of population dynamics, such as logistic growth for single species dy- 
namics in limited environment, Lotka-Volterra dynamics for predator-prey 
interactions, etc. 

• Self organized behavior. One expect to observe collective phenomena such 
as trophic chain formation or species formation. 

• Response to external stimuli, in particular to environmental changes such 
as cycle of seasons or human intervention. 

In classical population dynamics the building block are the species, and the 
interactions among them. Since we want to study the self-organization of an 
ecosystem (including species formation) we take as a building block the single 
individual. 

In our schematization the individual is identified by its genotype x, which is 
represented as a fixed length string of L bits: the genotype space is a Boolean 
hypercube of L dimensions, and mutations correspond to displacements in this 
space. On the other hand, a genotype identifies a strain of individuals. 

Individuals are able to survive according with a fitness function, which also 
takes into account the interactions with other individuals in the environment. 
Natural selection however does not act directly on the genotype, but rather 
on the resulting phenotype g(x), which can be considered a function of the 
genotype x* Generally, the phenotypic space is simpler than the genotypic 
one, according with the number of morphological characters considered. In our 
simple ecosystem model, g{x) is simply the fraction of ones in xJ In this case 
the phenotypic space is one-dimensional. The smoothness of the fitness function 
is related to the resolution required in genotypic space: if one clusters together 
a species into a single point then the fitness function can be quite rough. Since 
we are interested in the phenomenon of species formation, we require a smooth 
fitness landscape. 

*The assumption that the phenotype is a single-valued function of the genotype implies that 
we are not considering polymorphism (the fact that two cells with the same genotype can 
have different morphologies) nor age structure. 

^One can assume that in each locus there are two alleles of a given gene: for the "good" 
allele and 1 for the bad one 
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Finally, there is the real space. We shall describe in Section 2 a one- 
dimensional cellular automaton (CA) model. However, we shall consider only 
the limit of very long interaction length, i.e. global coupling. This simplification 
allows us to separate the complexity of the dynamics in genotypic space from 
spatial pattern formation. 

Similar systems have been introduced in order to investigate the phenomenon 
for which a phenotypically favored strain can loose its predominance due to a 



high mutation rate (error threshold). Bob In these works the population size 
is kept constant; recent preliminary results □ shows that if the population size 
is allowed to fluctuate (limited by an external constraint) another transition, 
called mutational meltdown Btj can be observed. In this case the whole popula- 
tion vanishes while not changing its distribution. No direct competition among 
individuals is considered. 

We are mainly interested in the problem of species formation due to inter- 
individual competition, in the limit of very small mutation rates. For this reason 
(and also due to the smoothness of the static fitness function), in our model 
the error threshold transition is not observable. Moreover, we do not impose 
any limit on the population size: individuals compete for free space and this 
automatically limits the size of population. The free space limitation translates 
into a logistic-like equation for the whole population size, and this furnishes a 
simple illustration of the mutational meltdown transition (see Section 3) . 

Analytic approximations can be obtained if one takes into consideration 
only the simpler phenotypic space, as reported in Appendix A. We are able to 
compute the speciation threshold, for which the population distribution splits 
into several separated peaks also for a very smooth fitness landscape and the 
mutational meltdown threshold. 

We developed an optimized computer algorithm for the original model, see 
Appendix B. The results of the simulations are reported in Section 3; the spe- 
ciation transition appears, for a choice of parameters consistently with our an- 
alytical results. In this more "realistic" case, the population distribution is not 
at all trivial, exhibiting coexistence of several quasi species at the same distance 
from the fittest strain. 



2. The cellular automaton model 

Let us consider an early ecosystem, populated by haploid- 1 " individuals. Each 
individual occupies a cell of a lattice in an one dimensional space; the size of 
the lattice is N sites. Each individual is identified by its genetic information 
(genotype), that we model as a base two number x of L bits. The distance 
in the genotypic space is defined in terms of the number of mutations needed 
to connect (along the shortest path) two individuals. We shall consider only 



* Single copy of genetic material, thus non sexually reproducing. 
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point mutations (0 <-> 1), occurring with probability \x independently of the bit 
position. Thus the genotypic distance d(x, y) between strains x and y is simply 
their Hamming distance (number of different bits). The mutation probability 
W(x,y) is 

W(x,y) = i/ (x ' v) 0--») L ~ d(x ' y \ 

which for vanishing mutation rate /U — > can be written in a quasi-diagonal 
form 

W(x,y) = n \id{x,y) = l 

W(x,x) = 1-Lfj, 

W(x, y) = otherwise. 

Given a genotype x, its phenotype is represented by a function g(x), which 
maps the genotypic space into the phenotypic one. In this paper we shall con- 
sider a very simple mapping, g(x) = d(x,0), i.e. the phenotype is proportional 
to the number of ones in the genotype. 

This automaton has a large number of states, one for each different genome 
plus a state (*) representing the empty cell. The evolution of the system is given 
by the application of two rules: the survival step, that includes the interactions 
among individuals, and the reproduction step. 

Survival: An individual Xi = x\ ^ * at time t and site i, i = 1, . . . , N, has 
a probability ir of surviving per unit of time. It is reasonable to assume this 
probability to depend only on phenotypic characters. The survival probability 
7r = n(H) is expressed as a sigma-shaped function of the fitness function H:§ 

e 0H 1 1 
n(H) = TT - m = - + - tanh(/3tf), (1) 

where (3 is a parameter that can be useful to modulate the effectiveness of 
selection. We always use (3=1. We define the fitness H of the strain Xi in the 
environment x = {x\, . . . , x* N } as 

1 N 

H(x it x) = h(g( Xi )) + J(9(xi),9{xj))- (2) 

i=i 

The fitness function is composed by two parts: the static fitness h(g(xi)), 
and the interaction term l/NY^j=iJ{g{ x -i)-,g( x j))- The matrix J define the 
chemistry of the world and is fixed; the field h represents the fixed or slowly 
changing environment. A strain x with static fitness h{g(x)) > represents 
individuals that can survive in isolation (say, after an inoculation into an empty 
substrate), while a strain with h(g(x)) < represents predators or parasites 

§Our choice of the fitness function does not consider the reproductive efficiency. 
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that requires the presence of some other individuals to survive. The interaction 
matrix J specifies the inputs for non autonomous strains. 

We assume that the static fitness h(u) is a linear decreasing function of u 
except in the vicinity of u — 0, where it has a quadratic maximum: 

h(u) = ho + b(l--- —V) ( 3 ) 
\ r 1 + u I r J 

so that close to u = one has h{u) ~ ho — bu 2 /r 2 and for u — ► oo, h(u) ~ 
ho + b(l — u/r). Thus, the master sequence (in Eigen's language) is located at 
x = 0. 

The matrix J mediates the interactions between two strains. For a classi- 
fication in terms of usual ecological interrelations, one has to consider together 
J{u,v) and J(v,u). One can have four cases: 

Jiu, v) < J(v,u) < competition 

J~(u,v) > J{v,u) < predation or parasitism 

J~(u,v) < J{v,u) > predation or parasitism 

J~(u,v) > J{v,u) > cooperation 

Since the individuals with similar phenotypes are those sharing the largest 
quantity of resources, the competition is stronger the more similar their pheno- 
types are (intraspecies competition). This implies that the interaction matrix 
J has negative components near the diagonal. We do not include here neither 
familiar structures nor sexual mating between genetically akin individuals nor 
other kind of competition or cooperation. 

We have chosen following form for the interaction matrix J: 

J{u,v) = -JK(^] (4) 



with the kernel K given by 



TV 



K(r) — cxp 

V a 

i.e. a symmetric decreasing function of r with K(0) = 1. The parameter J and 
a control the intensity and the steepness of the intraspecies competition, respec- 
tively. We shall use a Gaussian (a = 2) kernel, for the motivations illustrated 
in Appendix A. 

The survival phase is thus expressed as: 

• If Xi 7^ * then we get, with probability n(H(xi,x)), x\ = Xi, otherwise 

x[ = * 

• Else x\ = Xi = * 
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Reproduction: The reproductive phase can be implemented as a rule for 
empty cells: choose randomly one of the neighboring cells and copy its state; if 
it is different from the empty state then apply mutations by reversing the value 
of one bit with probability ji. 

One can notice that the effective reproduction rate does not only depend on 
the survival probability of the individual, but also on total availability of empty 
cells. 

One can have an insight of the features of the model by a simple mean field 
analysis. Let n{x) be the number of organisms with genetic code x, and rt* the 
number of empty sites, 

n* + n{x) = N. 

X 

We denote with m the relative abundance of non-empty sites: 
m = V] n ( x ) /N =1- n*/N. 

X 

The sums do not run over the empty cell state (x ^ *). We can express the 
fitness function H (and thus the survival probability 7r) in terms of the number 
of individuals n(x) in a given strain or in terms of the probability distribution 

p(x) = n(x)/mN 

H(x,n) = h(x) + J{ x ,v) n (y); ( 5 ) 

y 

H{x,p) = h(x) +m^2 J(x,y)p(y). (6) 
y 

The average evolution of the system will be governed by the following equa- 
tions, in which a tilde labels quantities after the survival step, and a prime after 
the reproduction step: 

n(x) = ir(x,n)n(x), (7) 

n'(x) = «W + ^E^>^)- ( 8 ) 
y 

Using the properties of W, 

y x 

and summing over x in Eqs. (0) and (]|), we obtain an equation for m: 

4 ^7r(x,n)n(x) = iriW (9) 

X 

y 



N 
N 
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i.e. 

m = rri7r(2 — mW), (11) 

where 

7f = — i- V 7r(s, n)n(x) = V 7r(£c, p)p(x) 
miV ^— ' 

a; a; 

is the average survival probability. 

The normalized evolution equation for p{x) is: 

Tr(x,p,m)p(x) + (1 -m7f)V V^(a;,y)7r(y,p,m)p(y) 
7r(2 — m7r) 



Notice that Eq. ( |ll| ) is a logistic equation with n as control parameter. The 
stationary condition, (m! = m), is 

2tF - 1 

m=^3-. (13) 

7T 

One observes extinction if 7f < 1/2. The decrease of 7f can arise from a 
variation of the environment (notably h(x)) or from an increase of the mutation 
rate /x, which broadens the distribution p(x). This last effect corresponds to 
the mutational meltdown, for which the population vanishes while continuing 
to exhibits a quasi-species distribution. Since the total population m multiplies 
the competition term in Eq. (fL2|), one cannot observe coexistence of species due 



to competition near the mutational meltdown transition. From Eq. (11) one 
could expect a periodic or chaotic behavior of the population; however, since if 
is always less than one, the asymptotic dynamics of the population m can only 
exhibit fixed points. 

We are mainly interested in the asymptotic behavior of the population in 
the limit fi — ► 0. Actually, the mutation mechanism is needed only to define 
the genetic distance and to populate an eventual niche. The results should not 
change qualitatively if more realistic mutation mechanisms are included. 

Let us first examine the behavior of Eq. ([jjj) in absence of competition 
(J = 0) for a smooth static landscape and a vanishing mutation rate. This 
corresponds to the Eigen model,aB with a smooth fitness landscape. Since it 
does not exhibit any phase transition, the asymptotic distribution is unique. 
The asymptotic distribution is given by one delta function peaked around the 
global maximum of the static landscape, or more delta functions (coexistence) if 
the global maxima are degenerate. The effect of a small but finite mutation rate 
is simply that of broadening the distribution from a delta peak to a bell-shaped 
curve 13 (quasi-species). 

While the degeneracy of maxima of the static fitness landscape is a very par- 
ticular condition, we shall show in Appendix A that in presence of competition 
this is a generic case. 
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Figure 1: Static fitness h, effective fitness H, and asymptotic, distribution p 
for the phenotypic model studied in Bagnoli and Bezzi (1997) cJ , analogous to 
Eq. (@), 



For illustration, we report in Figure [l] the asymptotic distribution of the 
population and the static and effective fitness for a similar model □ in which 
the genotypic space is approximated by a phenotypic one (see Appendix A for 
details). The effective fitness H is here almost degenerate, since /i is greater than 
zero and the competition effect extends on the neighborhood of the maxima, 
and this leads to the coexistence. 



3. Speciation and mutational meltdown in hypercubic genotypic space 

We have developed an optimized code (reported in Appendix B) for the sim- 
ulation of the original model. We use the following easter ^(/^representation for 
quasi-species in hypercubic space: starting from the origin of axis, we perform 
a step of a fixed length i?o with an angle nir/(2(L — 1)) — ir/8 if the nth bit 
(0 < n < L) of genotype x has value one. In this way one locates the master 
sequence (all zeros) at the origin; the strains with one bad gene, distributed ac- 
cording to the bad gene position at distance Rq\ the strains with two bad genes 
at an approximate distance of 2Rq, and so on. An example of the resulting 
hypercube for L = 4 is shown in Figure 0. 

We are interested in computing the critical values of parameters for the tran- 
sition between one single quasi-species to coexisting quasi-species (speciation). 



^We acknowledge D. Stauffer for having suggested this name. 
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Figure 2: The representation of the Boolean hypercube for L = 4 

We obtain from the approximate analysis of Appendix A that the crucial 
parameter in the limit /U — > is the quantity G = (J/R)/(b/r), which is the 
ratio of two quantities, one related to the strength of inter-species interactions 
(J/R) and the other to intra-species ones (b/r). We observe ,in good agreement 
with the analytical approximation see Eq. (|l5|), if Gm > G c (r/R) several quasi- 
species coexist, otherwise only the master sequence quasi-species survives. Here 
m is the asymptotic average population size that is usually close to one at the 
transition point. The approximate analysis also shows that this transition does 
not depend on the mutation rate in the first approximation. In Figure || a 
distribution with multiple quasi-species is shown. 

We can characterize the speciation transition by means of the entropy S of 
the asymptotic distribution, 



which increases in correspondence of the appearance of multiple quasi-species. 

In Figure || we characterize this transition as an increase of the entropy as 
function of Gm. We can locate the transition at a value Gm ~ 2.25, analytical 
approximation predicts G c (0.1) ~ 2.116. The entropy however is quite sensible 
to fluctuations of the master sequence quasispecies (which embraces the largest 
part of distribution), and it was necessary to average over several (15) runs in 
order to obtain a clear curve; for larger values of [i it was impossible to charac- 
terize this transition. A quantity which is much less sensitive of fluctuations is 




X 



10 F. Bagnoli & M. Bezzi 




.i'i- .* ~-' H 

« - 

.A 



Figure 3: Quasispecies in hypercubic space for L — 12. The smallest points 
represent placeholder of strains (whose population is less than 2 • 10~2), only 
the larger dots corresponds to effectively populated quasispecies; the size of the 
dot is proportional to the square root of population. Parameters: fi = 10~ 3 , 
ho = 2,b= 1CT 2 , R = 5, r = 0.5, J = 0.28, N = 10000, L = 12. 



the average square phenotypic distance from the master sequence g{x) 2 

SO) 2 = ^2,9(xfp{x). 



In Figure || (left) we characterize the speciation transition by means of g(x) 2 , 
and indeed a single run was sufficient, for [i = 10~ 3 . For much higher mutation 
rates (/i = 5 10~ 2 ) the transition is less clear, as shown in Figure || (right), 
but one can see that the transition point is substantially independent of /i, as 
predicted by the approximate theory, Eq. (|T^). 

Another interesting phenomenon is the meltdown transition, for which the 
mean field theory predicts extinction if w < 1/2, see Eq. (|T^). In Figure || 
we report the result of one simulation in which the extinction is induced by 
the increase of the mutation rate /i. One can notice that the transition is 
discontinuous, m jumps to zero from a value of about 0.15, and that the critical 
value of 7T is larger that the predicted one. This discrepancy can be caused by 
fluctuations, due to the finiteness of population. 

4. Conclusions 

We have studied a microscopic model of a simple ecosystem that exhibit the 
mutational meltdown effect and speciation phenomena. The size of the popu- 
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Figure 4: The speciation transition characterized by the entropy S as a function 
of the control parameter G m. Each point is an average over 15 runs. Same 
parameters as in Figure |[ varying J. 

lation is not hold constant; we found that in the mean filed approximation this 
quantity is ruled by a logistic equation, with the average fitness of population as 
control parameter. The model includes the competition among individuals, and 
this ingredient is considered fundamental for the speciation phenomenon in a 
smooth fitness landscape. This transition does not depend on the mutation rate 
provided that this rate is small. We are able to obtain analytical approximations 
for the onset of both transitions. 
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Appendix A: Analytical approximations 

Some analytical results can be obtained by considering the dynamics only in 
the phenotypic space. Let us consider the case of the phenotype that depends 
only on the number of bits (say, good genes) in the genotype, i.e. a highly degen- 
erate phenotypic space. One should introduce the multiplicity factor (binomial) 
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Figure 5: Independence of the speciation transition by the mutation rate. The 
transition is characterized by the average square phenotypic distance g{x) 2 of 
distribution p(x), as a function of the control parameter G m. Each point is a 
single run. Same parameters as in Figure ||[ varying J with fi = 10~ 3 (left) and 
[i = 5 1Q- 2 (right). 



of a given phenotype, which can be approximated to a Gaussian; however if 
one works in the neighborhood of the most common chemical composition, the 
multiplicity factors are nearly constants. Another reason for not using the mul- 
tiplicity factor is that we have not yet been able to derive analogous results with 
it. 

An instance of an application of a similar (sub-)space in the modeling of the 
evolution of real organisms is given by a repeated gene (say a tRN A gene) : a 
fraction of its copies can mutate, linearly varying the fitness of the individual 
with the "chemical composition" of the genets. This degenerate case has been 
widely studied (see for instance Alves and Fontanari (1996)liil); Another example 
is given by the level of catalytic activity of a protein. A non-degenerate space has 
also been used for modeling the evolution of RNA viruses on HeLa cultures.E3 

From now on we shall indicate with x both the phenotype and the genotype, 
and consider it as an integer number. To maintain a bit of the original multi- 
plicity, we extend the range of x to negative values, while keeping the master 
sequence at x — 0. 

We compute from Eq. (12) the values of parameters that allow the coexis- 
tence of different species. We look for a solution p(x) formed by the sum of 
delta peaks (fj, — > limit) centered at y^. 



P( x ) = ^2"fk5{x - y k ) = 

k k 
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Figure 6: Meltdown transition characterized by m and tt as a function of //. 
Here J = 0.3, h = 0.4, b = 0.35, r = 0.5, R = 5, N = 2000, L = 8. 



The weight of each quasi species is 7fc, i.e. 



J p k (x)dx = ik, 
The evolution equation for pk is: 



L-1 



fc=0 



Pk(x) = ^^Pk(x) 

TT 

The stability condition of the asymptotic distribution (p' k (x) = Pk(x)) is 
either Tr(yk) = tt = const (degeneracy of maxima) or pk{x) = (all other 
points). In other terms one can say that in a stable environment the fitness of 
all individuals is the same, independently on the species. 

The position yk and the weight jk of the quasi-species are given by Tr(yk) — 
tt = const and 3tt{x) / dx\ yk = 0, or, in terms of the fitness H, by 



L-1 



Vk ~ Vj 
R 



jj = const 



R 
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Let us compute the phase boundary for coexistence of three species for two 
kinds of kernels: the exponential one (a — 1) and a Gaussian one (a = 2). 
Numerical simulations show that the results are qualitatively independent on 
the exact form of the static fitness, providing it is a smooth decreasing function. 

Due to the symmetries of the problem, we have one quasi-species at x = 
and two symmetric quasi-species at x = ±y. Neglecting the mutual influence 
of the two marginal quasi-species, and considering that h'(0) = K'(0) = 0, 
K'(y/R) = —K'(—y/r), K(Q) = J and that the three-species threshold is given 
by 7o = 1 and 71 = 0, we have 

5(l-!) -K(z) = -1, i+K'(z)=Q, 
\ r I r 

where z — y/R, f = r/R and b = b/J. We introduce the parameter G — f/b — 
(J/R)/(b/r), that is the ratio of two quantities, one related to the strength of 
inter-species interactions (J/R) and the other to intra-species ones (b/r). In the 
following we shall drop the tildes for convenience. Thus 

r — z — mG exp \ ) = —mG, mGz 01-1 exp ( 

\ a I \ a 



Where m can be obtained from Eq. (|13|), 



2W — 1 e@ h ° 

~^T- with ?F = 7r (0) = TT -^- ) (14) 



and thus 



m = 1 — e 



-2f}h 



For a = 1 the coexistence condition never holds, except for G to = 1 and 
r = 0, i.e. a flat landscape (b = 0) with infinite range interaction (R = 00). 
Thus we suppose that the speciation transition is not present also for less steep 
potentials, such as power laws. 

For a = 2 the coexistence condition is given by 

z 2 — (mG + r)z + 1 = 0, mGz exp f — — ^ = 1. 

One can solve numerically this system and obtain the boundary G c (r) for the 
coexistence. In the limit r — > (static fitness almost flat) and (3ho ^> 1 (i.e. 
to ~ 1) one has 

G c (r) ~ G c (0) - r, (15) 



with G c (0) = 2.216 .... Thus for G > G c (r) we have coexistence of three or 
more quasi-species, while for G < G c (r) only the fittest one survives. The limit 
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(3h ^> 1 (to ~ 1) is not a restrictive condition from a theoretical point of view 
, in fact we can always stay in this approximation modulating (3; but we get it 
almost constant and equal to one, therefore there is a shortage of empty cells 
and the evolution can take much longer times. 

Appendix B: Monte Carlo Algorithm 

We describe here the essentials of the implementation of the model in FOR- 
TRAN (although we use C language). 

The implementation of the model can be done in a direct way, but since the 
coupling due to competition is global, the simulation time grows as TV 2 , where 
N is the number of individuals present in the environment. A way of speeding 
up a little the simulation is that of performing the computation of the fitness 
using Eq. (||) instead of Eq. (||) . 

The cellular automaton space is the vector integer env(0:N-l). The sur- 
viving individuals always occupy the first M positions (starting from 0): in- 
sertions always occurs at position M (M is thereafter incremented), and deleted 
individuals are overwritten by the genotype in position M-l (M is thereafter 
decremented) . 

We also use three other vectors, integer strain(0 : N-l) , integer distr(0:L2-l) 
and real fit(0:L2-l), where L is the genome length and L2=2**L. The first 
vector contains NS entries corresponding to each instance of a different genome 
in the environment. This vector is needed to perform sums over all genotypes 
without scanning env. The vector distr contains the number of instances of a 
given genome in the environment, and the vector fit its survival probability. 

The static fitness h and interaction matrix J are stored in the vectors real 
h(0:L2-l) and real J(0 : L2-1 ,0 : L2-1) , which are filled at the beginning. The 
central loop of the evolution algorithm is the following: 
C 

C assume that strain and distr are already OK 

C and compute fit 

C 

do i = 0, NS-1 
ig = strain(i) 
fit(ig) = 
do j = 0, NS-1 
jg = strain(j) 

fit(ig) = fit(ig) + J(ig, jg)*distr(jg) 
end do 

fit(ig) = fit(ig)/N + h(ig) 
fit(ig) = exp(f it(ig))/(l+exp(fit(ig))) 
end do 

C 
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C clear strain and distr 
C 

do i = 0, NS-1 

ig = strain(i) 

distr(ig) = 
end do 
NS = 

C 

C survival 
C 

i=0 

do while (i .It. M) 
r = rnd(iseed) 

if (fit(env(i)) .It. r) then ! don't survives 

env(i) = env(M-l) 
M = M-l 

else ! survives 

if (distr (env(i) ) .eq. 0) then ! first instance of genome 

strain (NS)=env(i) 

NS = NS+1 
end if 

distr (env(i) ) = distr (env(i) ) + 1 
i = i+1 
end if 
end do 

C 

C reproduction 
C 

Ml = M 

do i=Ml, N-l 

j = int (rnd(iseed) *N) 

if (j < Ml) 

if (rnd(iseed) .It. mu) then ! reproduction 

env(M) = ieor(env(j) ,2**(int(rnd(iseed)*L)) ! this is a XOR 

else 
env(M)=env(j ) 

end if 

if (distr (env(M) ) .eq. 0) then ! first instance of genome 

strain(NS)=env(M) 

NS = NS+1 
end if 

distr (env(M)) = distr (env(M) ) + 1 
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end if 
M=M+1 
end do 

The program has been implemented on a CRAY T3E and on a cluster of 
Linux machines using MPI. Since the code is not well parallelizable, due to the 
long range interactions and on the updating scheme, we have parallelized on the 
control parameters and on different runs. In other words we have launched a 
copy of the program in parallel on a different CPU, and the results have been 
collected using MPI. In this way also a cluster of machines with relatively slow 
connections (ethernet) can be used as a supercomputer. 
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